'''
Created on Mar 1, 2011

@author: oabalbin

The imput file is a genelist produced according to the vcfCodingsnps standar.
#name   chrom   strand  txStart txEnd   cdsStart        cdsEnd  exonCount       exonStarts      exonEnds        name    score   name2   cdsStartStat    cdsEndStat      exonFrames

BED format:
   1.  chrom - The name of the chromosome (e.g. chr3, chrY, chr2_random) or scaffold (e.g. scaffold10671).
   2. chromStart - The starting position of the feature in the chromosome or scaffold. The first base in a chromosome is numbered 0.
   3. chromEnd - The ending position of the feature in the chromosome or scaffold. The chromEnd base is not included in the display of the feature. For example, the first 100 bases of a chromosome are defined as chromStart=0, chromEnd=100, and span the bases numbered 0-99. 

The 9 additional optional BED fields are:

   4. name - Defines the name of the BED line. This label is displayed to the left of the BED line in the Genome Browser window when the track is open to full display mode or directly to the left of the item in pack mode.
   5. score - A score between 0 and 1000. If the track line useScore attribute is set to 1 for this annotation data set, the score value will determine the level of gray in which this feature is displayed (higher numbers = darker gray). This table shows the Genome Browser's translation of BED score values into shades of gray:
      shade                                                      
      score in range      166     167-277     278-388     389-499     500-611     612-722     723-833     834-944  945
   6. strand - Defines the strand - either '+' or '-'.
   7. thickStart - The starting position at which the feature is drawn thickly (for example, the start codon in gene displays).
   8. thickEnd - The ending position at which the feature is drawn thickly (for example, the stop codon in gene displays).
   9. itemRgb - An RGB value of the form R,G,B (e.g. 255,0,0). If the track line itemRgb attribute is set to "On", this RBG value will determine the display color of the data contained in this BED line. NOTE: It is recommended that a simple color scheme (eight colors or less) be used with this attribute to avoid overwhelming the color resources of the Genome Browser and your Internet browser.
  10. blockCount - The number of blocks (exons) in the BED line.
  11. blockSizes - A comma-separated list of the block sizes. The number of items in this list should correspond to blockCount.
  12. blockStarts - A 
'''

from collections import defaultdict

def read_known_genes_file(ifile,ofile,print_genes=True):
    '''
    This script does not preserve the isoform structure of different genes
    It only prints every unique gene or exon based on their location
    '''
    inputfile=open(ifile)
    outfile=open(ofile,'w')
    inputfile.next()
    gname, chr, strand, start, end, cdsstart, cdsend, exoncount, exonstart, exonends, name, score = \
    0,1,2,3,4,5,6,7,8,9,10,11
    padding=0
    
    genedict=defaultdict(list)
    exondict=defaultdict(list)
    for line in inputfile:
        fields=line.strip('\n').split('\t')
        if len(fields[chr].split('_')) > 1:
            print fields[chr]
            continue

        if print_genes:
            bedl=[fields[chr], fields[start], fields[end], fields[gname]+'|'+fields[name],fields[score],fields[strand]]
            loc=",".join([fields[chr], fields[start], fields[end]]).replace(',','&')
            if loc in genedict:
                continue
            else:
                outfile.write(",".join(bedl).replace(',','\t')+'\n')
            
            genedict[loc].append(fields[gname]) 
                        
        else:
            exons_starts=fields[exonstart].strip('').split(',')
            exons_ends=fields[exonends].strip('').split(',')
            
            for estart, eend in zip(exons_starts,exons_ends):
                eloc=",".join([fields[chr], estart, eend]).replace(',','&')
                if estart=='' and eend == '' or (eloc in exondict):
                    continue
                
                exondict[eloc].append(fields[gname])
                estart, eend = str(int(estart) - padding), str(int(eend) + padding)
                bedl=[fields[chr], estart, eend, fields[gname]+'|'+fields[name],fields[score],fields[strand]]
                outfile.write(",".join(bedl).replace(',','\t')+'\n')



def read_known_genes_file2(ifile,ofile,print_genes=True):
    '''
    It keeps all the isoforms for a particular gene
    It keeps all exons for the different gene isoforms
    '''
    
    inputfile=open(ifile)
    outfile=open(ofile,'w')
    inputfile.next()
    gname, chr, strand, start, end, cdsstart, cdsend, exoncount, exonstart, exonends, name, score = \
    0,1,2,3,4,5,6,7,8,9,10,11
    genedict=defaultdict(list)
    
    for line in inputfile:
        fields=line.strip('\n').split('\t')
        
        if print_genes:
            bedl=[fields[chr], fields[start], fields[end], fields[gname]+'|'+fields[name],fields[score],fields[strand]]
            outfile.write(",".join(bedl).replace(',','\t')+'\n')
                        
        else:
            exons_starts=fields[exonstart].strip('').split(',')
            exons_ends=fields[exonends].strip('').split(',')
            
            for estart, eend in zip(exons_starts,exons_ends):
                if estart=='' and eend == '':
                    continue
                
                bedl=[fields[chr], estart, eend, fields[gname]+'|'+fields[name],fields[score],fields[strand]]
                outfile.write(",".join(bedl).replace(',','\t')+'\n')

    
def sort_bed_file():
    '''
    It sorts the bed file according to chromosome and 
    in the same order as the ref genome of ucsc. chrM, chr1...,chr22, chrX, chrY
    '''
    pass

def build_exome_for_samtools(ifile,ofile):
    '''
    This script takes a bed file, e.g. exome
    and outputs the position file required for samtools 
    mpileup and bcftools
    chr1    66999824        67000051        NM_032291|SGIP1 0       +       1
    '''
    inputfile = open(ifile)
    outputfile= open(ofile,'w')
    i=0
    for l in inputfile:
        fields=l.strip('\n').split('\t')
        if fields[0]=="track db=\"hg19\" name=\"Human All Exon 50Mb targets\"":
            continue
        chr, start, end=fields[0],int(fields[1]),int(fields[2])
        while start <= end:
            outline=[chr, start]
            outputfile.write(",".join(map(str,outline)).replace(',','\t')+'\n')
            start+=1
    
            
    
    
'''
file_name='/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/CCDSGene.B37.txt'
ofile=('/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/CCDSGene.B37.bed')
efile=('/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/CCDSGene.B37.exons.bed')
'''


'''
file_name='/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/refGene.B37.txt'
#ofile=('/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/refGene.B37.bed')
#efile=('/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_30Mb_Kitv2/refGene.B37.exons.pad0.bed')

ofile=('/exds/users/oabalbin/projects/exomes/informative_genes_annotation.bed')
efile=('/exds/users/oabalbin/projects/exomes/informative_genes_annot_exons.bed')

print_genes=False
read_known_genes_file(file_name,efile, print_genes)
'''


# After this Intersect probes with exome. 
# Create a script with the bedtools line. Remember to count the probes. -a exome, -b probe_set

#ifile='/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/refGene.B37.exons.SureSelected.bed'
#ofile='/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/refGene.B37.exons.SureSelected.pos'
#ifile='/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_30Mb_Kitv2/refGene.B37.exons.SS30M.pad50.bed'
#ofile='/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_30Mb_Kitv2/refGene.B37.exons.SS30M.pad50.pos'

#ifile=('/exds/users/oabalbin/projects/exomes/informative_genes_annot_exons.bed')
#ofile=('/exds/users/oabalbin/projects/exomes/informative_genes_annot_exons.pos')


ifile='/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/SureSelect_All_Exon_50mb_with_annotation_hg19.bed'
ofile='/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/SureSelect_All_Exon_50mb_with_annotation_hg19.pos'
build_exome_for_samtools(ifile,ofile)
